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We introduce horizon pretracking as a method for analysing numerically generated spacetimes of 
merging black holes. Pretracking consists of following certain modified constant expansion surfaces 
during a simulation before a common apparent horizon has formed. The tracked surfaces exist at 
all times, and are defined so as to include the common apparent horizon if it exists. The method 
provides a way for finding this common apparent horizon in an efficient and reliable manner at the 
earliest possible time. We can distinguish inner and outer horizons by examining the distortion of 
the surface. Properties of the pretracking surface such as its expansion, location, shape, area, and 
angular momentum can also be used to predict when a common apparent horizon will appear, and 
its characteristics. The latter could also be used to feed back into the simulation by adapting e.g. 
boundary or gauge conditions even before the common apparent horizon has formed. 

PACS numbers: 04.25.Dm 04.70.Bw 



I. MOTIVATION 



In a spacetime that contains coalescing binary black 
holes, and given suitable gauge conditions, a common 
apparent horizon forms some time after the event hori- 
zons have merged. Similarly, in a spacetime containing 
a collapsing star, an apparent horizon forms some time 
after the event horizon has formed. In numerical calcu- 
lations it is important to locate apparent horizons. On 
one hand, one can extract interesting physical informa- 
tion from an apparent horizon, such as the object's mass 
and spin, for instance via the dynamical horizon formal- 
ism lHHD- On the other hand, certain numerical methods 
such as excision boundaries 0.0] use the apparent hori- 
zon to identify the region within the black hole. Know- 
ing the common horizon location also makes it possi- 
ble to adjust the dynamics of the spacetime via so-called 
"horizon-locking" gauges |5, 6], as well as gauges which 
attempt to control the location of the black hole 13]. At 
late times, the shape and oscillations of the apparent 
horizon is an effective indicator of the angular momen- 
tum of the final black hole, and has been related to its 
quasi-normal mode ringing |8]. 

From a practical standpoint, it is the apparent hori- 
zons which are more interesting than the event horizon 
at a given instant of time, since the latter is a globally de- 
fined quantity whose location can only be determined 
once the entire future development of the spacetime is 
known. The apparent horizon, on the other hand, is 
defined by the outermost marginally trapped surface 
and can be found locally in time. Further, it is always 
contained within the event horizon |9] when the strong 
energy condition holds, making it a "safe" estimator 
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of the black hole location for use by excision bound- 
aries 

A typical simulation scenario involves a pair of bi- 
nary black holes, identified by disconnected horizons in 
the initial data surface, evolving until they come close 
enough that a common horizon forms around them. The 
common apparent horizon appears instantaneously as 
a surface enveloping the two individual horizons. This 
is different from the situation for event horizons, where 
the disconnected lobes "grow together" to form a single 
connected surface |8, 9]. For a given dynamical evolu- 
tion, it is not known a priori when and where a com- 
mon apparent horizon first appears, as this depends on 
both the initial data as well as the particular slicing of 
the spacetime. 

Fast apparent horizon finders, which essentially solve 
a nonlinear elliptic equation, require a good initial guess 
for both the location and shape of the surface |12J. For 
horizons which are present in the spacetime from the 
initial time, it is usually sufficient to use the previously 
detected surface as an estimate for the horizon location 
at some small time interval later. The initial data con- 
struction itself often provides a good first guess for in- 
dividual horizons. The common horizon, however, ap- 
pears instantaneously at some late time and without a 
previous good guess for its location. In practice, an esti- 
mate of the surface location and shape can be put in by 
hand. The quality of this guess will determine the rate 
of convergence of the finder and, more seriously, also 
determines whether a horizon is found at all. Gauge 
effects in the strong field region can induce distortions 
that have a large influence on the shape of the common 
horizon, making them difficult to predict, particularly 
after a long evolution using dynamical coordinate con- 
ditions. As such, it can be a matter of some expensive 
trial and error to find the common apparent horizon at 
the earliest possible time. Further, if a common appar- 
ent horizon is not found, it is not clear whether this is 
because there is none, or whether there exists one which 
has only been missed due to unsuitable initial guesses — 
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for a fast apparent horizon finder, a good initial guess is 
crucial. 

In this paper, we present a reliable method for de- 
termining the first common horizon to appear in a dy- 
namical spacetime evolution. The method provides a 
close to optimal initial guess for the horizon surface and 
reliable predictions of its physical properties. Roughly 
speaking, horizon pretracking | 13] involves searching the 
spacetime for the surface of minimal constant expansion 
which envelops the source. For a spacetime without a 
common horizon, this surface will have positive expan- 
sion. Assuming that a common horizon does form after 
some time, the expansion of the pretracking surface will 
gradually decrease to zero with time, at which point the 
common horizon has formed, by definition. 

The rationale for this method is that it is much more 
reliable and efficient to "track" a surface than to find an 
entirely new surface, the difference being that a good 
initial guess is provided by the previous value of the 
tracked surface. In the past, a method such as pretrack- 
ing might have been too computationally expensive to 
be implemented as a practical solution. However, care- 
ful consideration of the problem of finding constant ex- 
pansion surfaces has resulted in algorithms which typi- 
cally solve the system in a matter of seconds Il l3.ll3tll4ll . 
Thus the required search for the minimal expansion sur- 
face becomes feasible. 

Pretracking has some similarities to minimisation fl5l 
Il6ll or curvature fP^l flow methods for apparent hori- 
zon finding. These methods solve a minimisation prob- 
lem or a parabolic equation to find horizons, starting 
from a large sphere, and can give a reliable answer as to 
whether a common horizon exists. Their disadvantage 
is that these methods are extremely slow and cannot be 
applied at every time step. Pretracking is faster since 
it "tracks" instead of minimising or "flowing" at each 
time. Additionally, the pretracking surfaces can provide 
valuable information about the current state of the sim- 
ulation. 

In Section UTI we briefly discuss surfaces of constant 
expansion. Section ITTT1 defines the notion of a pretrack- 
ing surface, and outlines a number of algorithmic details 
involving the parametrisation of the surface, and the bi- 
nary search method for the minimal constant expansion 
surfaces. Finally in Section HVl we present applications 
of the technique to binary black hole mergers. 



II. APPARENT HORIZONS AND CONSTANT 
EXPANSION SURFACES 

An apparent horizon (AH) is the outermost 
marginally outer-trapped surface in a spacelike hy- 
persurface I. The AH satisfies the conditions 0/« = 
and 0( n ) < 0, where 0^) and 0/„) are the expansions 
of the outgoing and ingoing null normals of the surface. 



These can be calculated from the ADM variables via 

= +V i s i + K ij {s 1 si -yi) and (1) 

(n) = -Vij + Kijijsi-yi). (2) 

Here yu and Ku are the three-metric and extrinsic curva- 
ture of the spacelike hypersurface, and s, is the spacelike 
normal to the surface within the spacelike hypersurface. 

We describe the surface explicitly through a func- 
tion h(9,(p) which specifies the radius of the surface 
as function of the angular coordinates 9 and (p about 
some origin. The spacelike normal is then given by 
S; = VjF/l VF|, where F(r, 6, cp) = r - h(9, <p) with the 
coordinate radius r, hence F(r,9,(p) = indicates the 
surface location. This explicit representation is conve- 
nient but not necessary for surfaces with S 2 topology, 
and is not restrictive in practice H3l . 

The notion of pretracking requires specifying a gen- 
eral family of surfaces for spacetimes which may not 
have an apparent horizon, but which contain the AH as 
a surface if one exists. For instance, a generalisation to 
the concept of apparent horizons are the constant expan- 
sion (CE) surfaces, which are defined by the condition 

e (/) = c, (3) 

with C a constant over the surface Q . The apparent 
horizon is a CE surface with C = 0. For Ku — 0, 
the CE surfaces are identical to constant mean curvature 
(CMC) surfaces, which are defined by V,s ! = C. In 
an asymptotically flat spacetime with zero extrinsic cur- 
vature, the CMC surfaces foliate the hypersurface in a 
neighbourhood of infinity |18]. Thus a parametrised 
family of CMC surfaces can be specified at large radii, 
though they may not exist in the strong field interior 
of a spacelike slice. In numerical experiments with iso- 
lated sources, we have found that CE surfaces form 
a similarly useful foliation at large distances even for 
Kij 7^00. 

A simple algorithm for pretracking can be described 
as follows: Define a parametrised family of CE surfaces 
hp (9, (p) which envelop the source, via 

h v {9,<$>) : © {e) [h p ] =C P , (4) 

where p labels each member. Here we introduce the no- 
tation [h] to indicate that a quantity is a functional of the 
surface function h(9,4>). The surface with the smallest 
positive expansion Cp forms the "best guess" for the lo- 
cation of the AH. For time evolutions of vacuum space- 
times, this surface is expected to change continuously, 
and can be followed until the expansion of the surface 
reaches zero, at which point it represents the earliest 
common AH. 

A difficulty arises from the above definition of the 
pretracking surface hp which introduces some practi- 
cal complications. While the expansion is zero at the 
location of the AH, it also decreases to zero asymp- 
totically towards spacelike infinity, while being posi- 
tive in between. This is illustrated in Fig. for flat 
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Figure 1: Expansion and the quantity r© vs. coordinate 
radius r for flat space in Minkowski coordinates and for a 
black hole in Brill-Lindquist and Kerr-Schild coordinates. Both 
= and r© = define the location of the apparent hori- 
zon, which is at different radii in different coordinate systems. 
While reaches a maximum at a finite radius, r© increases 
monotonically in spherical symmetry. 



space in Minkowski coordinates and for a black hole in 
Brill-Lindquist and Kerr-Schild coordinates. For hori- 
zon finding we are more interested in the surfaces most 
closely enveloping the source. Thus, the initial guess 
must be chosen inside the maximum of if the method 
is to avoid converging to one of the asymptotic sur- 
faces. To remove this ambiguity, a more suitable func- 
tion would increase monotonically outside the common 
AH. In practice, one therefore has to use a modified min- 
imisation criterion, as described below. 



III. HORIZON PRETRACKING 
A. Method 

The pretracking method attempts to find a surface 
corresponding to a "generalised apparent horizon". 



That is, we would like to find a horizon-like surface 
which envelops a compact source, even at times when 
an actual common horizon does not exist, but which 
reduces to the apparent horizon when it is present. 
The generalised surface is found among a specified 
parametrised family, {hp}. The pretracking method cor- 
responds to finding the "smallest" member of the fam- 
ily, i.e., the value of the parameter p such that the surface 
hp (9, cp) satisfies 



H[hp] = const, 
G[hp] = min, 



(5) 
(6) 



where the shape function H[h] is a generalisation of the 
expansion 0^y and the goal function G[h] specifies what 
we mean by "smallest". Note that H maps to a function, 
while G maps to a scalar. 

The parametrised family {hp} should be defined in 
such a way as to contain the common AH if it exists 
in the given slice. The shape function, H, specifies the 
surfaces that we want to find. We allow here for a 
more generic set of surfaces than only CE surfaces. The 
goal function, G, defines a suitable notion of "closeness" 
among the surfaces H p to an actual apparent horizon. 
As such, it is useful to define it in such a way that it 
evaluates to zero for a surface which is a common AH, 
negative within the AH, and monotonically increasing 
outside the AH. 

We generalise the concept of CE surfaces for a number 
of pragmatic reasons. For instance, very distorted sur- 
faces are difficult or impossible to represent accurately 
in terms of a function h{9, (p), so we have found more 
rounded surfaces work consistently better. We have ex- 
perimented with the following shape functions: 



Htlh] 
H r [h] 
H r2 [h] 



® {e) [h], 

{i) [h]-h, 
{e) [h]-h 2 . 



(7) 
(8) 
(9) 



As above, h(9, <p) is the shape of the surface, i.e., the sur- 
face's coordinate radius at the coordinate angles 9, (p. 
The first definition, H\, is simply the surface's expan- 
sion, so that the condition H\ [h] = const specifies CE 
surfaces. We find empirically that the shape functions 
H r and H r z lead to surfaces with a less distorted coordi- 
nate shape. (This can be seen e.g. in the lower left hand 
graph of Figure[6|) The shapes H r and H r 2 are in general 
not CE surfaces, but they tend to CE surfaces for large 
radii and close to the apparent horizon. 

Note that H r and H r 2 are not defined in a covariant 
manner, because their definition depends on the coordi- 
nate shape. While we would prefer covariant shapes, we 
find that the listed non-covariant shapes work much bet- 
ter in practice and are computationally more efficient. 
Empirically, especially the shape function H r leads to a 
reliable pretracking method. All of the above surface 
shapes are apparent horizons if and only if the shape 
function is zero everywhere. 



4 



The goal function G [h] must have a minimum for the 
surface hp that we want to call "closest" to being an AH. 
That is, the goal function has a minimum which is lo- 
cated at the AH when one exists. We considered the fol- 
lowing goal functions: 



G H [h] := H[h], (10) 
G Hr [h] :=H\h]-h, (11) 
G r [h) := h, (12) 

where the over-bar 7 denotes the average over the sur- 
face. We take this average as an average over all grid 
points. One could alternatively define a covariant aver- 
age that takes the two-metric into account, but we did 
not consider the advantage of this to be worth the addi- 
tional c ompl exity in the equations for the Jacobians (cf . 
Section Till Bl and Appendix lA) . For a solution h, we re- 
quire not only a specific value of G, but also that H[h] be 
constant over the surface, so that e.g. Ghj- = implies 
H = everywhere. 

The combination of the shape function Hi and the 
goal function Gfj leads to the simple algorithm men- 
tioned in the previous section. As already discussed, 
the difficulty with this choice is that even in spherical 
symmetry, the goal function is not monotonic, but rather 
reaches a maximum at a finite radius. Thus, this method 
is not expected to converge unless an appropriate ini- 
tial surface, within the maximum, is chosen. One way 
to resolve this ambiguity is to redefine the shape func- 
tion, for instance via JSJ or @. For either of these surface 
functions, the goal function Gh increases monotonically 
(at least in spherical symmetry) and is zero on the AH. 
The lower graph of Fig. shows the goal function Gh 
for the shape function H r , i.e., the quantity rO. 

The combination of the shape function H\ and the 
goal function Gur also ensures a monotonically increas- 
ing goal function in spherical symmetry, with the ad- 
ditional advantage that all pretracking surfaces are CE 
surfaces, and thus covariantly defined. Unfortunately 
we find empirically that this combination does not work 
as well in p ractice as the combination of H r and Gh (see 
Section^). 

The goal function G r measures the average coordinate 
radius of a given surface, so that the method selects the 
surface with smallest average coordinate radius satis- 
fying H[h] = 0. Note however, that the AH does not 
usually form at the innermost surface defined by H (as 
shown in in Section ITW Thus this goal function cannot 
reliably be used for pretracking, although it can still be 
very useful to study the behaviour of the innermost CE 
surface. 

Finally, we note that the goal functions Gj-jr and G r 
are not covariantly defined, as they depend on the coor- 
dinate radius of the surface. This is not an issue, as the 
value of the goal function itself has no relevance other 
than the value G [h] = 0, which by definition is the co- 
variantly defined AH. 



B. Surface Finding 

Finding a surface h(8,ep) that satisfies one of the 
above conditions H[h] = const under a constraint 
G[h] = p for a given value p is comparable to finding 
an apparent horizon. Efficient implementations of ap- 
parent horizon finders have been developed in recent 
years I13tll4ll . and it is a straightforward extension to al- 
low them to look for alternate similarly defined surfaces, 
such as those defined by the H* listed above. Current 
efficient AH finders interpret the surface defining equa- 
tion as a nonlinear elliptic equation; they use a Newton- 
Raphson iteration method to linearise and then a Krylov 
subspace method to solve it. 

The overall elliptic equation that defines the pretrack- 
ing surface is 



H[h](0,4))-H[h] + G[h]-p = 0, (13) 

which is fulfilled if and only if H [h] is constant over the 
surface and if G[h] = p. This is because the first term 
H[h] is the only term that varies over the surface, and 
therefore it has to be constant for the equation to be ful- 
filled. In this case, the expression H[h] — H[h] vanishes, 
which then implies G[h] = p. Discretised, this equation 
becomes 

H i( h j) ~lL H k(hj) + G(hj) - p = 0, (14) 

JV k 

where i, j, k label grid points, and N is the total number 
of grid points on the horizon surface. The Jacobian of 
this equation is 

JV k 

where /,•; = dHj/dhj is the Jacobian of the surface shape 
function H,, and djG = dG/dhj. For example, when H 
is the shape function H\, the Jacobian ]u is the Jacobian 
of the original apparent horizon equation. 

There is one crucial problem that appears when one 
attempts to solve equation 1141 . The Newton-Raphson 
method requires the Jacobian of the equation to be 
sparsely populated, otherwise the method will be pro- 
hibitively expensive. 1 Unfortunately, the Jacobian \15\ 
is densely populated due to the term Y,k Jkj> which is 
nonzero for all values of the indices i and j. For N grid 
points, this increases the storage requirements by a fac- 
tor 0(N), and turns the 0(IV 3 / 2 ) runtime cost of the 
solver into a cost of 0(N 3 ). This is in general not ac- 
ceptable. 



1 A non-local discretisation of the horizon shape, e.g. an expansion in 
spherical harmonics, would not require a sparse Jacobian. 
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We arrive at a sparse matrix by extending the vector 
space of the solution /z, by one element. Instead of ilii , 
we write the equivalent system 



Hi(hj) 
\-G(hj) 



C 

- V 



o, 



(16) 



where we have introduced one additional unknown C 
and one additional equation that determines it. The Ja- 
cobian of this equation is 



Ji 



N 



djG 



+1 



(17) 



This Jacobian is now sparsely populated. In addition to 
the already sparsely populated Jy, it has one additional 
fully populated row and column, so that it has 0(N) 
nonzero entries out of 0(N 2 ) total. The system M6t can 
therefore be solved efficiently. 

We list the full expressions for the Jacobians for the 
different shape and goal functions in AppendixlAl 



finding an apparent horizon with the advantage of a 
good initial guess provided by the previous search. This 
takes usually less than a second per step on current 
hardware, independent of the resolution of the under- 
lying spacetime. Various trade-offs can be made to im- 
prove efficiency. For instance, pretracking need not be 
carried out at every time step of an evolution, but rather 
at longer intervals. We usually determine p* to only a 
moderate accuracy at each pretracking time. 

In order to find the common apparent horizon as early 
as possible, we follow up with a regular apparent hori- 
zon search using the pretracking surface as initial guess, 
which is a very good initial guess near the time when 
the common AH forms. This finds the common appar- 
ent horizon even when pretracking is not accurate, and 
also allows us to pretrack in larger time intervals. The 
inaccuracy in the horizon guess provided by pretrack- 
ing comes mainly from the fact that failure to find a pre- 
tracking surface for a certain parameter value p does not 
necessarily indicate that the surface does not exist. Since 
we usually compromise some accuracy for efficiency in 
the pretracking search, it may be that we did not iterate 
long enough or did not start with a good enough initial 
guess for this value of p. 



C. Pretracking Search 

At each time step during an evolution, pretracking 
consists of determining the parameter p which selects 
a member h of the family of CE surfaces {hp} that min- 
imises the goal function G[h] . A convenient initial guess 
for p is either the value from the last pretracking time or 
a manually specified guess, for instance corresponding 
to a large sphere. Because the equation defining the sur- 
face is nonlinear, it is also necessary to specify an initial 
guess for the surface shape h itself, again either from a 
previously determined surface or manual specification. 

We have modified an apparent horizon finder so that 
we can specify the desired shape function H, the goal 
function G, the desired parameter value p, and an initial 
guess ho for the surface. The result is either a surface 
h that satisfies H[h] — const and G[h] = p, or a flag 
indicating that no such surface could be found. 

The pretracking surface generally exists only for a cer- 
tain parameter range p > p* , i.e., only for values of p 
above some critical parameter p*. During the search, 
we keep track of an interval [p m i n , p m ax] that indicates 
which values of p are known to fail and succeed, respec- 
tively. We start with a trial value of p that was the result 
of the last pretracking time, and depending on whether 
the corresponding surface can be found or not, we set 
either p m ; n or p max - We then increase or decrease p in 
large steps until we find the other end of the interval. Fi- 
nally, we use a binary search within the interval to find 
p* to a given accuracy Ap, and call the resulting surface 
the pretracking surface at this time. 

Each step of the above algorithm, i.e., each check of 
a parameter value p, is computationally equivalent to 



D. Apparent Horizon Tracking 

When a common apparent horizon first forms, we 
find that it bifurcates into an inner and an outer hori- 
zon. This has, for example, been described by Thorn- 
burg 1 14] (cf . his Figures 3 and 4). If the apparent hori- 
zon world tube is spacelike, then the spacelike hyper- 
surface L can intersect ("weave") the world tube in ar- 
bitrary ways, and generically, such intersections will oc- 
cur in pairs. This can be seen e.g. in |19, Fig. 3], which 
outlines the trapped region in a spherical collapse. At 
early times, the trapped region's boundary is spacelike. 
The outer boundary, which corresponds to the apparent 
horizon, approaches the event horizon at late times. The 
structure of individual, inner, and outer horizons is also 
visible in |20, Fig. 1] or |21, Fig. 1]. 

We demonstrate this for the Brill-Lindquist data in 
the next section. For black hole evolutions, one usually 
chooses gauge conditions that have the effect of mak- 
ing the inner horizon quickly shrink in coordinate space, 
while the outer horizon initially grows and then stays 
roughly constant in radius. 

Of course, we would like to ensure that horizon track- 
ing follows the outer horizon, as it is the one of physical 
interest to observers in the exterior spacetime. However, 
unless special measures are taken, the horizon finder 
will select a branch at random. We observe that shortly 
after bifurcation the outer horizon is generally less dis- 
torted in coordinate shape in a binary black hole merger. 
Before the inner and outer horizons have clearly sep- 
arated in coordinate space, the tracker can also jump 
between the two branches. We can make the horizon 
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tracker select the outer horizon by smoothing its initial 
guess. Similarly, we can make it select the inner hori- 
zon through a more distorted initial guess. Addition- 
ally by slightly enlarging or shrinking the initial guess, 
the outer or inner horizons are preferred, respectively. 
We implement this by modifying the initial guess, fen,, 
before each horizon search by the prescription 



fen 



[l-f)-g-h+f-g-h , 



(18) 



where 7 is the average over the surface. The smoothing 
factor, f, should be positive for smoothing and negative 
for roughening. We find that / = 0.25 is a good smooth- 
ing choice for the particular data we examined. The en- 
largement factor g must be positive, and will grow (or 
shrink) the surface if it is larger (or smaller) than one. 
We find that g — 1 . 1 is usually a reasonable value to en- 
sure that an outer horizon is found. On the other hand, 
/ = —0.2 and g = 0.8 selects the inner horizon in our 
example of a head-on collision discussed in the next sec- 
tion. 



IV. APPLICATIONS 

We applied the algorithms described in the previous 
sections to 3D simulations of (1) an axi-symmetric head- 
on collision and (2) coalescing binary black holes with 
angular momentum. For the time evolution we used 
the conformal-traceless (BSSN) formalism introduced 
in 1 22, 23, 24]. The particular implementation we use, 
including gauges (1 + log lapse and T-driver shift) and 
boundary conditions, is described in [25]. The binary 
black hole initial data are calculated via the puncture 
method [26]. Spins of the surfaces were measured via 
the dynamical horizon formalism, as described in @]. 
Our time evolution code is implemented in the Cactus 
framework |22]. 

Pretracking was performed using extensions to 
Thornburg's AHFinderDirect horizon finder 11411 . A 
number of different shape and goal functions were 
tested. These are listed in Table Qj along with corre- 
sponding labels which will be referenced in the plots in 
this section. 



A. Head-on Collision 

As a first test case, we studied a simple axi- 
symmetric head-on collision of two black holes with 
time-symmetric Brill-Lindquist initial data 1 28], evolved 
in 3D. The black holes have a mass parameter m = 1, 
and are located on the z-axis coordinate locations Z = 
±1. The total ADM mass of the system is therefore 
M — 2. The resolution for the particular runs performed 
here is Ax = 1 /8, and the outer boundaries were placed 
at coordinate x max = y max = 4 and z max = 5. We chose 
a pretracking surface and apparent horizon resolution 
of 5° and a pretracking accuracy of Ap = 10~ 4 . 



Label 


Shape Goal tracked 
function function horizon 


e 


Hi 


G H 


outer 


r© 


H r 


G H 


outer 


r 2 


H r2 


G H 


outer 


r-bar 


Hi 


G[ir 


outer 


smallest r 


Hi 


G r 


outer 


inner r& 


H r 


G H 


inner 



Table I: Algorithms described in this paper are listed according 
to their use of shape functions and goal functions. The meth- 
ods tested in this paper correspond to the various shape and 
goal functions defined in Section llll Al The names listed in the 
leftmost column are used to label curves in the figures in this 
section. 



As the singularities are punctures, we evolved them 
without excision by ensuring that their locations are 
staggered between grid points |25]. For the lapse, a, we 
used a 1 + log slicing condition, starting from a = 1, 
and used normal coordinates, ft' = throughout the 
evolution. The simulation was performed with an ex- 
plicit octant symmetry. The model is particularly useful 
as a test case, as it begins with disconnected apparent 
horizons, but very soon (after t fa 1) forms a common 
horizon, requiring minimal resources in terms of com- 
putation time or memory for the given resolutions and 
boundary locations. 

Figure |2 shows various pretracking results from this 
evolution. The top left graph shows the average of the 
expansion versus time for pretracking surfaces from 
the various algorithms listed in Table 0] The algorithms 
0, r0, r 2 &, and r-bar all find the common apparent 
horizon at identical times, around t = 1 . (However, the 
algorithm r-bar unfortunately fails in more realistic 
simulations; see e.g. Figure [6]) The algorithm smallest r 
fails to find the common horizon. It locates the CE sur- 
face with smallest average radius, a surface which has 
a positive average expansion even after a common ap- 
parent horizon has formed. This surface actually inter- 
sects the pretracking surfaces defined by the other algo- 
rithms, and seems "smallest" due to the averaging of the 
radius over the surface which includes a narrow throat. 

The top right plot shows the average of the expansion 
of the ingoing null normal ®( n ) which must be negative 
for apparent horizons. The lower two graphs show the 
average coordinate radius (left) and areal radius (right) 
versus time of the pretracking surface. After an initial 
transient, all of the pretracking algorithms behave sim- 
ilarly, with the exception of the smallest r, which fails as 
a pretracking surface. The gradual growth in coordinate 
radius of the common apparent horizon is caused by 
the shift condition; a zero shift is unsuitable for a long- 
term evolution since the black hole eventually encom- 
passes the grid. Note that the inner r0 algorithm tracks 
as desired the inner apparent horizon (as described in 
the previous section), which has a smaller radius. 
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Figure 2: Averages of the expansion ®, inner expansion &i n y coordinate radius r, and areal radius R for a head-on collision vs. 
coordinate time t. Except for the smallest r case, the common apparent horizon is found when = 0. Shown are the results for 
five different pretracking algorithms. The methods labelled O, r®, r 2 ®, and r-bar O track minimum surfaces of the corresponding 
quantity, as described in the main text. The smallest r algorithm tracks the CE surface with the smallest coordinate radius, and 
fails to converge to the apparent horizon. Finally, inner r® uses the same pretracking algorithm as r®, but locks onto inner instead 
of outer horizons. The roughly parabolically shaped region between the inner and outer horizons in the lower left hand graph 
contains trapped surfaces. 



Before a common horizon is formed, the different 
pretracking algorithms actually determine different sur- 
faces. Figure [3] shows the shapes of the pretracking sur- 
faces at t — 0,1, and 1.5 in one quadrant of the xz- 
plane. Although the three successful pretracking algo- 
rithms determine rather different surfaces initially, they 
have converged by t = 1, shortly before the first com- 
mon apparent horizon is found. At t = 1.5 the first three 
algorithms continue to track the same outer horizon sur- 
face, while the difference between the outer and inner 
apparent horizon (as tracked by the inner r& method) 
is clearly visible. The algorithm smallest r, which fails, 
tracks a very different surface. This surface is neverthe- 
less interesting because it is, up to numerical errors, the 
smallest CE surface that exists within the given slices. 



1. Evolution ofCE Surfaces 



It is interesting to study the behaviour of various CE 
surfaces independent of the pretracking algorithm. Fig- 
ure 01 compares how CE surfaces with constant areal 
radius R and CE surfaces with constant expansion 
evolve in time. Inside an inner cutoff, the CE surface 
equation has no solution. We are using normal coordi- 
nates, and constant areal radius surfaces grow in coor- 
dinate space. New common CE surfaces with smaller 
areal radii and smaller expansions begin to exist near 
the inner cutoff. One of these new surfaces is the com- 
mon apparent horizon with expansion = 0. 
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Figure 3: Shapes of the pretracking surfaces for various pre- 
tracking methods at different times. The coordinate time t = 
shows the initial data. The time t = 1 is shortly before the 
common apparent horizon is detected. Here the different pre- 
tracking algorithms already produce very similar results. At 
t = 1.5 a common apparent horizon has formed, and the inner 
and outer horizons have separated. The individual apparent 
horizons are also shown for comparison with the pretracking 
surfaces. The region between the inner and outer horizons in 
the lowermost graph contains trapped surfaces. 



Grid spacing 
Ax 


Parameter Areal radius 
p* R 


1/4 


0.05498 


3.880 


1/5 


0.05665 


3.893 


1/8 


0.05540 


3.887 


1/10 


0.05483 


3.889 


1/16 


0.05507 


3.885 


1/20 


0.05508 


3.882 



Table II: Convergence test for different spatial resolutions. The 
values seem to approach a certain limit when the resolution is 
increased, indicating that the limit limA.v^o V* does exist. 



2. Convergence Properties of the Method 

We used the above setup to test convergence of the 
pretracking code on the initial data set for the method 
rO. There are several numerical steps involved in pre- 
tracking, namely first the time evolution of the space- 
time, then an interpolation to the pretracking surface, 
then the solution of an elliptic equation on the pretrack- 
ing surface, and finally the pretracking search for the 
critical parameter value p*. Each of these steps may or 
may not converge to a certain appropriate order. The in- 
tention here is to test the accuracy of the pretracking it- 
erations, assuming that the remaining steps are already 
convergent, in order to demonstrate that this pretrack- 
ing search is a well-defined problem, i.e., that such a 
critical value p* does indeed exist independent of the 
resolution given appropriately defined shape and goal 
functions. 

The convergence of the interpolation and the appar- 
ent horizon finder has already been reported in I14f], and 
our modifications to make it find pretracking surfaces 
instead of apparent horizons do not affect this. We omit 
studying a time evolution by considering only Brill- 
Lindquist initial data. We only change the resolution 
Ax of the data given on the 3D hypersurface, and cal- 
culate the critical parameter p* as a function of Ax. We 
leave the resolution of the pretracking surface constant 
at 2.25°, and set the pretracking tolerance to Ap = 10~ 8 . 

The resolutions Ax, resulting critical parameters p*, 
and areal radii R, of the numerically determined sur- 
faces are given in Table |H1 The results converge for 
higher resolutions, but a fourth order convergence is not 
evident. This is likely caused by the fact that finding a 
surface depends not only on the desired goal function 
value p, but also on the initial guess Jiq for the surface. 
This makes it in practice very difficult to determine the 
critical parameter value p* rigorously. The given val- 
ues p* and R have therefore errors that are unrelated to 
the grid spacing, which makes a convergence test prob- 
lematic. We take here a pragmatic approach and remark 
that pretracking does in practice lead to a reliable detec- 
tion of the common apparent horizon immediately after 
it forms. The precise value of p* is here not of interest. 
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Figure 4: Left (top and bottom): Time evolution of various CE surfaces with given areal radii R. The graphs show CE surfaces 
with areal radii from 3.25 to 6 (from bottom to top). The upper left graph shows that these surfaces grow with time, the lower 
left graph shows that their expansion decreases with time. The largest CE surface is not found after t w 3.5 any more because it 
reached the outer boundary. The three smallest CE surfaces only begin to exist after some time. 

Right (top and bottom): CE surfaces with given expansions 0. The coordinate radii of the surfaces increase with time, and new 
surfaces come into existence near the inner boundary where the CE surfaces cease to exist. The common AH is one of these 
surfaces. For the smallest new CE surfaces also the corresponding inner surfaces (with the same expansion) are shown. 



3. Pretracking Efficiency 

We show in Fig. [5] the cost of pretracking. As de- 
scribed in Section ftll CI pretracking iterates over several 
surfaces until a surface with a certain minimum prop- 
erty has been found. Each such pretracking iteration 
is about as expensive as finding an apparent horizon. 
A "typical" simulation of our group runs in parallel on 
about 30 processors and needs about 20 seconds for each 
evolution time step. Finding an apparent horizon with 
a good initial guess takes about one second, i.e., about 
5% of the evolution step. The time spent in the appar- 
ent horizon finder is approximately independent of the 
resolution of the 3D simulation and the number of pro- 
cessors. The efficiency of our fast horizon finder is de- 
scribed in detail in 1 14]. 

In this simulation, the successful pretracking methods 



take about eight iterations during pretracking, and one 
iteration while tracking the horizon later. We find that 
it is in general sufficient to pretrack only every 10 iter- 
ations or fewer, and the cost of pretracking is therefore 
acceptable. 



B. Coalescing Black Holes with Angular Momentum 

As a further application for pretracking we consider 
an inspir ailing black hole setup. We emphasise that pre- 
tracking is of high importance in these cases, where the 
two black holes can stay apart for times of more than 
100 M (see e.g. ] 29]). Pretracking provides valuable in- 
formation about impending apparent horizon formation 
as well as the necessary initial guess for the apparent 
horizon finder, or it can (with some confidence) exclude 
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Figure 5: Number of pretracking iterations, which has the 
same cost as the same number of AH searches, for different 
pretracking methods at different times. The pretracking accu- 
racy was Ap = 10~ 4 ; this means that the value p* of the pre- 
tracking surfaces at each time is only accurate up to this error 
(see section lirm . 



that a common AH exists. 

The model we consider is the innermost stable circu- 
lar orbit as predicted in |30], which applies the effective 
potential techniques of |31] to puncture initial data |26]. 
Previous simulations with this dataset suggest that they 
perform around a half orbit before a common horizon 
first appears H^HHl- The initial data contain two punc- 
tures with a proper horizon separation of L = 4.99 M, 
angular momentum / = 0.78 M 2 and angular velocity 
fl = 0.17 /M, where M is the ADM mass of the sys- 
tem. This dataset was studied in detail in 1 33] . This 
model is the first in the sequence (QC-0) studied with 
the Lazarus perturbative matching technique in 1331 . 
The QC-0 model is attractive as a test-bed for pretrack- 
ing since a common apparent horizon is found at around 
t = 16.44 M, while the simulation continues accurately 
beyond t = 40 M. 

In this model, an initial guess for the common appar- 
ent horizon with a spherical form will find the common 
AH at the same iteration as pretracking if the chosen ra- 
dius is accurate to about 10%. For black hole simula- 
tions that start out with a larger separation, the initial 
guess usually has to be ellipsoidal in order to find the 
common AH at all. In that case, all axis lengths have to 
be guessed correctly to the above accuray. 

The evolutions were performed in bitant symmetry, 
i.e., with a reflection symmetry about z = plane. We 
used a grid containing 256 x 256 x 128 points and a res- 
olution of Ax = 0.06 M, which places the outer bound- 
ary at only about 7.68 M for this test run. The evolu- 
tions were stopped at t = 22 M since we were not inter- 
ested in further tracking the evolution of the common 
apparent horizon and its ring-down. During the evolu- 
tion we used a co-rotation shift to keep the black holes 



from moving across the grid urn, which is essential 
for long-term stable evolutions. We used a 1 + log slic- 
ing and a f-driver shift condition |25|]. 

Figure [6] shows the evolution of the pretracking sur- 
faces found by different algorithms. The upper left 
plot shows the various pretracking surfaces at the ini- 
tial time. 2 It is obvious that the different algorithms are 
locating very different surfaces. Note that the covariant 
method r-bar has the most distorted shape. The bump 
in the waist of the r 2 surface at t = 5 M and t = 10 M 
is caused by insufficient resolution in the surface; our 
explicit h(8,(p) representation is not well adapted to a 
peanut-shaped surface. The surface evolution can still 
be tracked, however, and later approaches the correct 
AH shape. 

After t — 12.75 M, the method r-bar fails to locate a 
surface. It succeeds again at t = 16.79 M, but it has then 
locked onto a surface that is inside the AH. We assume 
that this is due to the larger coordinate distortion of this 
surface, which has a much more pronounced waist than 
the other methods. The individual horizons are lost af- 
ter t = 14.34 M. At t — 16 .44 M a common appar- 
ent horizon is found for the first time. The lower right 
hand graph compares different surfaces at that time and 
demonstrates that the succeeding pretracking methods 
have converged to the same surface. 

Figure |7| shows the average expansions and the areal 
radii of the pretracking surfaces. The time at which the 
common AH is found is marked with a vertical line. 
The different pretracking methods track different sur- 
faces during most of the evolution, but these surfaces 
converge about 1 .5 M before the common AH is found. 
The pretracking surfaces approach the common appar- 
ent horizon smoothly, and they accurately predict both 
its shape and its area. Pretracking gives here an early in- 
dication of the time at which the common AH forms. Af- 
ter the common AH has formed, the pretracking meth- 
ods r0 and r 2 lock onto the AH. The method r-bar 
0, which failed to locate a surface for some time, after- 
wards tracks a different surface. Both methods r0 and 
r 2 are thus reliable pretracking methods for this evolu- 
tion. 

We measured the angular momentum of the pre- 
tracking surfaces using the dynamical horizon frame- 
work |1]; our implementation is described in |2]. We 
show the results in Fig. [5] Note that this angular mo- 
mentum measure is a quasi-local quantity and therefore 
is defined on any surface, not only apparent horizons. 
However, our current implementation is restricted to 
cases that are close to having an axial Killing vector. In 
dynamical spacetimes this condition will only be satis- 
fied approximately. Yet in the case of this merger simu- 
lation both successful pretracking methods r0 and r 2 



2 The simplest method, O, did not find an appropriate pretracking 
surface in this case, and is not shown. 
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Figure 6: Time evolution of the different pretracking surfaces in the xy-plane for the inspiralling QC-0 model. Shown are 
individual and common apparent horizons and the pretracking surfaces found using different pretracking algorithms. A common 
horizon first appears at t = 16.44M, shown in the lower right plot. It does not show the r-bar surface because it could not be 
found at that time. 



predict the angular momentum of the common apparent 
horizon. 



We conclude that pretracking works as an effective 
analysis tool for this binary black hole collision. The 
time of merger is predicted by the rate in which the 
expansion of the pretracking surfaces approaches zero, 
and the pretracking surfaces give a good estimate for 
the shape, area, and angular momentum for the future 
common apparent horizon. 



V. CONCLUSION 

Horizon pretracking is an effective method for pre- 
dicting the location and shape of an emerging common 
horizon, as well as yielding important physical informa- 
tion during the course of a spacetime evolution. The 
simple cases tested here indicate that it is applicable to a 
variety of scenarios, such as black hole mergers, where 
the formation of a horizon is expected, but its location 
and shape are not known a priori. 

Various practical issues need to be taken into consid- 
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Figure 7: Average expansion and areal radius R for the in- 
spiralling black hole evolution. The average expansion for r 2 
is larger than can be shown on this scale. This is due to numer- 
ical errors, as explained in the main text. The method r-bar © 
fails after t = 12.75 M; it succeeds again after the common AH 
has formed, but tracks a surface that is not the AH. 



eration when choosing an appropriate surface family for 
pretracking. A family of CE surfaces parametrised by 
their expansion suffers some drawbacks, since the ex- 
pansion does not increase monotonically. However, a 
simple modification such as the multiplication by a co- 
ordinate radius leads to an effective practical algorithm. 
We have shown that several methods for defining the 
pretracking surface all lead to procedures which accu- 
rately determine the common AH at the time of its first 
formation. 

As the lifetimes of black hole evolutions with dynamic 
gauges are extended, this ability to predict the common 
horizon shape will become particularly important. Ap- 
parent horizon finders require a reasonable initial guess 
in order to converge. After a long period of evolution of 
separated black holes, the initial common horizon shape 
is difficult to predict. As a result, a common AH may not 
be found, or at least not at the earliest possible moment. 
This can lead to the impression that one is evolving 



„ 1-0 

M 

-> 0.8 
E 

I 0.6 

CD 



% 0.2 
0.0 



r© 
r 2 
r-bar 
x common AH 



E 0-4 S ^ ++4 * t +++++ +++++ *++ + Jp<. 



xxx 




10 
time [M] 



15 



20 



Figure 8: Angular momentum / of the pretracking surfaces 
versus time, measured with the dynamical horizon frame- 
work. The angular momentum of the common AH increases 
directly after it has formed, indicating that it accretes both en- 
ergy and angular momentum initially. In the few cases where 
a value for J could be found for the r-bar pretracking sur- 
faces, they agree well with the r& method. In some cases, no 
angular momentum could be measured because no approxi- 
mate Killing vector could be found on the pretracking surface 
at that time. 



separated black holes when they have actually already 
merged. Valuable physical information, which can be 
determined from the shape and oscillations of the com- 
mon horizon, will be lost. On a practical level, finding 
the earliest common AH gives important information as 
to the causality of the spacetime and allows more effec- 
tive application of certain numerical techniques, such as 
singularity excision. 

We observe that new trapped surfaces form in pairs of 
an outer and an inner surface. A given initial guess for 
the horizon surface will find only one of these surfaces, 
which tend to diverge quickly in coordinate space. Since 
initially they are identical, however, an AH finding al- 
gorithm will converge to either the inner or the outer 
depending on its initial guess. This is particularly true 
if the common AH is found very early (as is the inten- 
tion of the pretracking method). We have given a simple 
algorithm that helps selecting the outer, generally less 
distorted surface, which is the astrophysically interest- 
ing apparent horizon. This algorithm does not depend 
on pretracking and is useful in general. 

We emphasise here that valuable information can be 
gained from the study of foliations of CE (and related) 
surfaces in a slice, beyond the commonly studied appar- 
ent horizons. The dynamical horizon formalism defines 
a quasi-local measure of the spin that can be evaluated 
on any closed surface, and the three-covariantly defined 
CE surfaces form a natural choice. Gauge conditions can 
also be designed to reduce distortions in such surfaces, 
or hold them in place on the grid, potentially simplify- 
ing the dynamics. 
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The pretracking method is efficient enough to be ap- 
plied as a regular analysis tool during large simulations. 
It is already being used in a variety of simulations of bi- 
nary puncture and thin sandwich data evolutions, and is 
finding further application in hydrodynamical collapse 
scenarios. 
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Appendix A: EXPRESSIONS FOR THE JACOBIANS 

We list here for the reader's convenience the equa- 
tions and Jacobians for the pretracking methods that we 
use in this paper. See table U for the combinations of 
shape and goal functions and the names that we use for 
these combinations. 

We assume that a desired goal function value p is 
given, and a corresponding surface h is to be determined 
by the combination of the shape and the goal function. 
For all methods, we list first the system of equations that 
defines the surface, and then give the Jacobian of this 
system. The intermediate quantity C is also determined 
through these equations, but its value is irrelevant. 

We denote the expansion of the surface h with 0,-, and 
the expansion's Jacobian as ]u = d@j/dhj. These expres- 
sions are e.g. given in 11311 or 11411 . is the Kronecker 
delta symbol. 

Method 0: shape function Hi, goal function Gh- 



Oj-C = 
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Method r 2 0: shape function H r %, goal function G# : 
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Method r-bar 0: shape function Hi, goal function Gh/- 
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Method smallest r: shape function H\, goal function 
G,: 

0,-C = (A9) 
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Method inner r@: same as r& above 
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